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Motivated by the recent rapid development of the field of quantum gases in optical lattices, we 
present a comprehensive study of the spectrum of ultracold atoms in a one-dimensional optical 
lattice subjected to a periodic lattice modulation. Using the time-dependent density matrix renor- 
malization group method, we study the dynamical response due to lattice modulations in different 
quantum phases of the system with varying density. For the Mott insulating state, we identify 
several excitation processes, which provide important information about the density profile of the 
gases. For the superfluid, the dynamical response can be well described in a local density approxi- 
mation. This simplification can be valuable in understanding the strong-correlated superfluid in a 
slow-varying harmonic potential. All these spectroscopic features of an inhomogeneous system can 
be used as a test for the validity of the Bose-Hubbard model in a parabolic trapping potential. 

PACS numbers: 03.75.Lm, 03.75.Kk, 05.30.Jp, 73.43.Nq 
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I. INTRODUCTION 

Following the rapid development of experimental tech- 
niques for the manipulation and detection of dilute ultra- 
cold atom g wide range of fundamental quantum 
many-body phenomena have been observed. Specifically, 
due to techniques including Feshbach resonances [lj and 
optical lattices @, bosonic systems loaded into a peri- 
odic lattice described by the Bose-Hubbard modelfj-ly] 
have been experimentally accessible both in the weakly 
and strongly interacting regimes with highly control- 
lable parameters, allowing for example the observation of 
the superfluid to Mott-insulator phase transition driven 
by quantum fluctuations Q- This achievement has pro- 
vided a new platform to study quantum many-body 
physics by virtue of the high degree of control and tuning 
available Q. 

Particularly rich quantum physics is to be expected 
in the context of quantum many-body physics far from 
equilibrium, but this remains largely unexplored at the 
moment. One main experimental difficulty lies in the lim- 
ited number of techniques of measurement in the strongly 
correlated system. In this paper, we will focus on one 
particular technique and point the way to the extraction 
of additional theoretical information from the raw data, 
the periodic lattice modulation approach by Stoferle et 
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al. @ They developed this technique to study the excita- 
tion spectrum of the bosonic system in an optical lattice. 
It acts at a probe with a specific frequency on the ul- 
tracold bosons and can be used to reveal the excitation 
spectrum of the system. More recently, this technique 
has been widely used in the dynamical contr ol [Tol a nd 
the realization of the quantum phase transition nf 13 [ in 
optical lattices. 

Previous theoretical and numerical studies on this ex- 
periment have shown how to extract important informa- 
tion about the system. Theoretically, the technique has 
been studied via perturbative methods in two limits, by 
a linear response analysis [id], [l5[ in the Mott insulating 
phase and by solvin g th e Gross-Pitaevskii equation in 
the superfluid regime [16(. A drawback of these perturba- 
tive methods is that they cannot be used to deal with 
the whole interaction regime. Numerically, the time- 
dependent density-matrix renormalization group tech- 
nique (t-DMRG) has been applied to simulate the ex- 
perimental setup in a quasi-exact fashion (l7l - [l9| . Basic 
features seen in the experiment could be reproduced suc- 
cessfully. 

Although these studies have opened a window on the 
understanding of the experimental observation 9] , further 
theoretical and numerical questions about the excitation 
spectroscopy arise mainly for three reasons. First, ultra- 
cold experiments are carried out in a harmonic trapping 
potential. Although it had been pointed out that this will 
induce spectral broadening and a shifting of peaks, a de- 
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tailed and quantitative study of this issue is still lacking. 
This is important to assess whether the presence of the 
trap qualitatively or quantitatively changes the behavior 
of the homogeneous system. Secondly, a direct compari- 
son to the experimental observation, which averages over 
many ID systems with different particle numbers and 
densities per ID tube, is not satisfactory unless the de- 
pendence on the density is taken into account. Therefore, 
building upon previous studies, this paper will generally 
focus on extracting new information from the spectrum 
due to the harmonic confinement and different densities. 



II. THE BOSE-HUBBARD MODEL 



Ultracold bosons in an optical lattice can be described 
by the Bose-Hubbard model@, HJ 



H = -J s T(b}b 



j 



(i) 

where bj and fij are the annihilation and number opera- 
tors on site j, respectively. The first term describes the 
hopping process between nearest neighbours, while the 
second one depicts the on-site interaction. The last term 
models the harmonic trapping potential Vj — Vt(j — jo) 2 
with Vt the curvature and jo the center of the system. 
In this paper, we are only interested in the absorption 
spectrum of ultracold bosons in a ID optical lattice. We 
assume the ID optical lattices all to be directed along the 
x-direction. We denote by V x and V± the laser strength 
along the x-direction and yz-directions respectively. For 
deep lattices, the hopping matrix element J and on-site 
interaction U can be approximated as [171. |20| 



J_ 




(2) 



III. METHODS 

A. Time-dependent perturbation (t-perturbation) 

If the modulation amplitude 5V is small and the sys- 
tem is in a deep Mott-insulating state (J<Ct/), the sys- 
tem can be understood within the framework of time- 
dependent perturbation theory where the Hamiltonian 
reads 



H[J(t),U{t)]^H + H'{t), 



(4) 



where H Q = H{t = 0). (For brevity U = U(t = 0) 
and Jo = J(t = 0) are used hereafter.) To keep the 
only term contributing to the excitations, we make a 
transformation 141 



(5) 



By further neglecting the time-independent term, we 
have 



L-l 



H'(t) = -Fj sinwt 53(StS i+1 + H.c), 



with the coupling constant 

/<91nJ d\nU 

Fj 



dV x 



dV x 



JnSV. 



(6) 



(7) 



V X =V 



This couplingconstant has been shown to be valid in the 
large U- limit [15[- 

In standard time-dependent perturbation theory, the 
transition probability is given by 



W mn (t) 
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(8) 



and 



U i — a, 



! v x vl 



(3) 



where E r is the recoil energy, a s is the s-wave scattering 
length, and A is the wavelength of the laser forming the 
optical lattice. The parameters used in the calculations 
are a s = 5.45 nm, A = 825 nm, V± = 3QE r . In order 
to investigate the absorption spectrum, one applies a si- 
nusoidal modulation of the x-direction laser strength V x 
starting at t = with frequency ui and amplitude 5V ', i.e., 
V x (t = 0) = Vb and V x {t > 0) = V + SVsinut^, and 
measures the energy absorbed by the system. The ab- 
sorbed energy is strongly frequency-dependent and gives 
information about the excitation spectrum. 



Here H' mri is the matrix element of H' between two eigen- 
states to and n of the unperturbed Hamiltonian Hq. So 
the energy absorbed is 

AE(t) = J2 h ^noW m0 {t). (9) 



B. t-DMRG 

Numerically, we use the t-DMRG method to study 
the time evolution of the system [22j-|24|. This method 
is a quasi-exact algorithm which allows for simulating 
real time evolutions of ID quantum many-body systems, 
which operates on a class of matrix product states [25T4271 ] . 
To begin with, a conventional finite-system DMRG algo- 
rithm is used to determine the ground state, \ip(t = 0)), 
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of the Hamiltonian at time t = 0, H(t — 0) for a system 
with L sites and N bosons. Then a full time evolution 
of the quantum state, is calculated with the t- 

DMRG algorithm based on a Trotter decomposition of 
time steps. We keep up to 200 states in the reduced 
Hilbert space in the algorithm. In order to reduce the 
error from the Trotter decomposition, we use a linear fit 
to extrapolate the results to Trotter time steps 5t — > 0. 
Convergence in the number of states kept has also been 
checked, such that on the time scales simulated the re- 
sults are quasi-exact. 

The full time dependence of the total energy reads 

E( t ) = m)\H(t)\m) 

^ (m\H \m) + (m\H'(t)\m)- m 

The main contribution to the energy transfer to the sys- 
tem comes from the first term as the time average of 
H'(t) oc smut vanishes. Thus, to get the absorption 
spectrum, we calculate the energy absorbed up to a given 
time t m , 

AE = (iP(t = t m )\H{t = 0)|V>(t = tm)) - E Q , (11) 

where Eq is the ground state energy This method is 
essentially equivalent to fitting the full time dependence 
of the energy curvefl7l|. 



IV. RESULTS 



(a) 




FIG. 1. (Color online) (a) Total energy absorbed as a function 
of the modulation frequency uj subjected to a small modula- 
tion 5V = O.OlVb- Both the results of t-DMRG and time- 
dependent perturbation theory are shown. The inset shows 
the initial density distribution of the system, (b) and (c) are 
schematic descriptions of the excitations due to the modula- 
tion. 



A. Absorption spectroscopy with Mott domains 

The Mott phase in the Bose-Hubbard model is char- 
acterized by short-ranged exponentially decaying corre- 
lations and commensurate filling. In a trapped system, 
there is no homogeneous Mott phase because of the har- 
monic trap. But for large enough U and suitable particle 
filling, the system may still have one or several Mott do- 
mains separated by superfluid domains in space fvA fl8j. 
The number of Mott domains and the filling in those do- 
mains depend on the average density of the system. In 
the following, we will study the absorption spectroscopy 
in the presence of Mott domains with varying densities. 

For the low density case, we consider a system with 
12 bosons in a deep optical lattice with Vq = 15E r 
and 36 sites. The curvature of the trapping potential 
is V t = 0.0123EV w O.OWo- The densit y profile of the 
system is shown in the inset of Fig. l(a)| where there is 
only one Mott domain with unit filling in the center of 
the system. Then we consider the absorption spectrum 
which is measured at time t m = 100h/E r and depicted 



in Fig. 1(a) In a homogeneous Mott phase, the absorp- 
tion spectrum is highlighted by a sharp peak at energy 
U 17]. This resonance, corresponding to a particle-hole 
excitation, is sharp since the excited states are almost 
degenerate. However, the situation is different after a 



trapping potential is applied. The trap introduces a dif- 
ference in potential energy between two neighboring sites 



AV(j,j±l)=V t [l±2(j-j )]. 



(12) 



As a rough estimate, the particle-hole excitation energy 
from site j to j + 1, as shown in Fig. |(lb)| will deviate 
from U by V(j, j + 1). The width of the peak is expected 
to be determined by the potential difference at the edges 
of the domain where the potential energy is maximized. 
According to the density profile shown in the inset of 



Fig. 1(a) the edge is at site 13 and site 24 where the 



potential difference is 0.17Uq. So the estimated width of 
the 1U peak is 0.34[/q, which coincides very well with our 



numerical results in Fig. 1(a) 



Another new feature in the spectrum is the small peak 
at haj s» 0.21C/o which corresponds to the particle-hole ex- 
citation where a particle hops to an empty site as shown 



in Fig. 1(c) This stems mainly from hopping from site 



13 to site 12 and site 24 to site 25, where the potential 
difference is 0.204[/o, in excellent agreement with numer- 
ics. 

As the particle number N increases, there will be sev- 
eral Mott domains with different fillings as shown in the 
inset of Fig. 2(a) where the particle number is 36 and 



the other parameters are the same as in the previous 
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(a) 




1 1.5 2 2.5 3 3.5 

Modulation frequency Rfl)/U 




FIG. 2. (Color online) Absorption spectrum (a) in a confined 
system with 36 particles. The inset of (a) shows the initial 
density profile, (b) and (c) are schematic descriptions of the 
excitations between two different Mott domains, (b) describes 
the second predominant absorption peak at tko = 1.75/7o in 
(a), while (c) is used to explain the small peak at ftu = 0.25[/o- 
From the inset one can calculate the chemical potential dif- 
ference at the Mott domain walls via Eq. I12I which gives 
AV = AV(J,j + 1) = AVQ.2, 11) = AV(25, 26) = 0.25(7 . 



case. Besides the broadened 1U peak, there are two more 
peaks: one centered at 0.25Uo, and another centered at 
1.75f/o- The former corresponds to the hopping of elec- 
trons between two Mott domains with different doping at 
the domain boundaries as depicted in Fig. [5J and its fre- 
quency is determined by the potential energy difference 
AV(12, 11) = AV(25, 26) = 0.25<7 - 

At a first glance, the peak at 1.75£/q originates from 
the hopping of bosons from the unit filled Mott domain 
into the doubly occupied region. This is qualitatively 
but not quantitatively accurate. If there is no exter- 
nal potential, the excitation energy is exactly 2U, in- 
dicating the non-unit filling of the system 17]. However, 
due to the parabolic potential, one has to take into ac- 
count the difference of the potential energy of the two 
site involved in the hopping process AV as shown in 
Fig. ~ 
shou 

and Fig. 2(c) involve exactly the same two sites in the 
system, so we have AV — 0.25£/o which is the peak en- 
ergy analyzed above. Then the correct position should 
be at 1.75 /To which is exactly the results in our numer- 
ical calculations. This shift has also been observed by 
Stoferle et al., who reported that there was a peak at 



2(a) Thus, the position of the excitation peak 



d be hu = 2U — AV. Both excitations in Fig. 2(b) 



about 1.9(7 in the Mott phase @. 

An interesting consequence of these observations is 
that by combining the broadening effect of the 1/7 peak 
and the shift of the 2U peak, one can directly determine 
U. In the case considered here, we find 



U 



(13) 



with Wi the width of the 1U peak and U2 the position 
of the 2U peak. This formula works because the broad- 
ening and shifting effect are caused by almost the same 
chemical potential difference AV. This is useful in the 
calibration of the on-site interaction parameter U. 

(a) 




1 1.5 2 2.5 

Modulation frequency nco/U 

(b) 



3.5 




0.5 1 1.5 2 2.5 3 3.5 

Modulation frequency nto/U 



FIG. 3. (Color online) For a weak modulation of a large 
system with 58 particles we have (a) the total energy absorbed 
vs. the modulation frequency. On the other hand, in (b) we 
show the spectroscopy of the system in Fig[2] under a strong 
modulation 8V = 0.2Vb- 



In order to further corroborate the connection between 
the energy shift and the density profile, we performed an 
additional calculation in a Mott system involving triple 
occupancies shown in Fig. 3(a) In this high-density sys- 
tem, there exist particle-hole excitations at the various 
Mott domain walls. Since the chemical potential differ- 
ence is spatially varying, the shifts away from 2U are 
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different. Therefore, we can see the AV peak and 2U 
peak both split into two peaks due to this difference. 

To compare our results with experiments where the 
modulation amplitude is up to 20% of the lattice depth, 
we also carried out a simulation for a large modulation 
(Fig. |3(b)[ ). Here the breakdown of perturbation the- 
ory indicates the saturation effect in a real system. Due 
to this effect and the relatively large trap curvature the 
splitting of the 1U peak is pronounced. We can still iden- 
tify the positions of the excitations from the shifts. For 
example, the peak at 1.17?7o is related to the hopping 
process from site 14 to 13, or from site 23 to 24. Another 
important finding is that the positions of the AV peak 
and 2U peak are robust. Also, the saturation effect for 
these excitations is less significant than that of the 1/7 
resonance. 

Based on the findings above, the position of the 2U 
peak reveals important information of the density profile 
of the Mott system. For one thing, the number of the 
peaks indicates the number of Mott layers in the "wed- 
ding cake" structure. For another, the shift away from 
2U provides important information about the positions 
of the Mott domain walls. 



B. Absorption spectrum in a superfluid 

In this section, we turn to the superfluid state, where 
the interpretation of the excitation spectroscopy is less 
straightforward than in the Mott regime. Without loss 
of generality, we choose Vo — ^E r and SV — 0.2Vo in all 
the calculations in this subsection, leading to U/J s» 5. 
In contrast with the Mott insulating state where the hop- 
ping is substantially suppressed, both the parameters J 
and U play an important role in determining the main 
properties of the superfluid. What makes the situation 
more complicated is the external harmonic trap, which 
introduces inhomogeneity into the system. A simplifica- 
tion occurs nevertheless as we will show that one can map 
the absorption spectrum in a trapped system to the ho- 
mogeneous one by using the local density approximation 
(LDA). Mathematically, this means 

A£ trappcd H = / p(r)Ae homo (w, p{r)) dr, (14) 

where p(r) is the spatially dependent density, 
A£ tra PP cd (u) is the energy transferred as a function of 
frequency u in a trapped system, and Ae homo (w, p(r)) is 
the energy absorbed density in a homogeneous system 
with particle density p(r). 

To show this approximation really holds, we compare 
the exact spectrum of a trapped system with a result 
from the LDA (Fig. |4(b)[ )[2g. To simplify the calcula- 
tion of spectra wc first make an approximation on the 
density profile (Fig. |4(a)j ). Although this approximation 
seems rough, the resulting spectrum is in good agreement 
with the exact one. Then the approximate results can 
be calculated as A£ LDA =EiPW ehom °( w -P( ! )). where 



e homo (w, p(i)) is calculated in a homogeneous system with 
density p{i). 

From the comparison, it is clear that the LDA works 
extremely well despite the approximate density profile. 
From a physical point of view, the validity of LDA stems 
from the slow- varying density profile in the harmonically 
confined superfluid. Also the chosen interaction is away 
from the location of the phase transition. Thus, here the 
main effect of the parabolic trap is no more than intro- 
ducing a slow- varying inhomogeneity. On the other hand, 
LDA must fail in the Mott insulating phase since there 
exists non-trivial excitations at the boundaries between 
Mott domains, as we saw in the previous section. 



(a) 




o 1 1 1 1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 
Modulation frequency nco/U 



FIG. 4. (Color online) (a) Density profiles for the system. 
Both the exact (circles) and the approximate (squares) results 
are shown, (b) The absorption spectroscopy in the form of 
the exact t-DMRG results (circles) and LDA (squares). The 
lines are drawn to guide the eye. 



The significance of the LDA is that one can under- 
stand the properties of a trapped superfluid with the 
help of a homogeneous one. To further test the LDA 
and to study how the density influences the spectroscopy, 
we performed calculations for different particle numbers 
both in homogeneous and confined systems (see Fig. [SJ. 
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(a) Homogeneous 




1 1.5 2 2.5 3 3.5 
Modulation frequency nto/U 

(b)Confincd 
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1.5 2 2.5 3 3.5 
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FIG. 5. (Color online) Spectra with different particle number 
N in (a) a homogeneous system and (b) trapped one. The 
system size L = 36 is fixed. 



spectrum: in the low-density regime, it exhibits a two- 
peak structure; as the particle number increases, the two 
peaks merge and become a broad one. This subtle change 
is also an indication for the density of the quantum gases. 
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For the unconfined atoms, there is a sharp excitation 
between 2>U and AU for intermediate densities. This reso- 
nance can be interpreted as two particle-hole excitations 
with both particles on the same site [3]. Here, we find 
that this peak is very sensitive to the density of the quan- 
tum gas. As shown in Fig. 5(a)| when the density of the 
system is below 1/2, this excitation is negligible because 
it is rare for three particles to be nearest neighbours. As 
the density increases to unity, i.e., N — 36, the spectrum 
shows a sharp peak around 3.8J7. Strictly speaking, its 
excitation energy involves two parts. The first part is the 
excitation energy to the 3[/-Hubbard band, while the sec- 
ond one comes from the change of the kinetic energy for 
these two atoms from the delocalized ground state to the 
localized excited state. In a system with unit filling, the 
energy difference per atom between derealization and lo- 
calization is about 2 J. Therefore, the total energy gained 
for this type of excitations should be 3C7 + 2 x 2 J. In our 
case where J 0.2[/, the estimate 3.8U coincides with 
the t-DMRG result. 

When we consider the system in a harmonic trap (see 
Fig. |5(b)| , we find a similar density dependence for the 



FIG. 6. (Color online) The centers of the peaks in the spectra 
of Fig. [5] for (a) a homogeneous system and (b) trapped one. 
The values are obtained by fitting the curves with bimodal 
gaussian distributions. 

To demonstrate the resemblance of the two cases (ab- 
sence or presence of a trap), we also plot the centers of 
two peaks in the spectra as a function of the number of 
particles in Fig. [51 It is clear that the two peak centers 
shift towards each other as the particle number increases 
in both cases. For the first peak, it moves towards higher 
energy because hoppings between sites with differing fill- 
ings become important with more particles. For the sec- 
ond one, its excitation energy of the second peak involves 
3U for the Hubbard-band and 4 J for the kinetic energy. 
It shifts towards lower energy because the ground state 
for atoms is not completed delocalized for a high density, 
which compromises the energy gained from localization 
in the excited state. From this argument, it is expected 
that the peak would be located around 3C7, which is con- 
sistent with our numerical result. With further increased 
density, the 3U resonance and the excitations around 2U 
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merge together and form a broad continuum eventually. 

Thus, the density dependence of the absorption spec- 
trum in either a homogeneous or trapped system is basi- 
cally the same. To understand this in the framework of 
the LDA, a trapped system can be divided into several 
homogeneous region with slow- varying density. The en- 
ergy transferred A£' trapped is mainly determined by the 
high-density region since it maintains the most important 
weight in the integral of Eq. ((14)) . Then it is reasonable 
to map the trapped system to a homogeneous one. 

V. DISCUSSION AND CONCLUSION 

How can these observations on the absorption spec- 
troscopy in a trapped system be turned into tools for 
experimental analysis? It turns out that they lead to 
schemes to test the validity of the Bose- Hubbard Hamil- 
tonian and to calibrate the interaction parameters. 

First of all, our findings in the Mott regime are useful 
in determining the density profile of the cold gases in an 
optical lattice. As has been mentioned in Sec. IIV Al the 
broadening of the 1U excitation and the shift of the 2U 
resonance are both directly related to the density dis- 
tribution of the system. For example, the number of 
the 2U peaks indicates the number of the Mott domains, 
while their shifts away from 2U can help to calculate the 
boundaries of the domains. Also, it is straightforward to 
test the Bose-Hubbard model by comparing the results 
from the spectroscopy with those from the time-of-flight 
imaging technique. Another important application in the 
Mott phase is calibrating the parameter U, which has 
proven to be difficult in a deep lattice 29]. The simplest 
way to do this with the spectroscopy is to fit the first 
strong 1U peak. To improve the accuracy, one can take 
into account the 2U resonance to recalibrate. For in- 
stance, the U is directly related to the width of the first 
peak and the position of the second peak according to 
Eq. (fKT)) . The advantage of this method is that no other 
fit parameters are needed. 

For the superfluid case, the main concern is how to 
make use of the LDA. Since one can map the absorption 
spectrum of a confined system to a homogeneous one, it 



suggests that the basic features of the superfluid, includ- 
ing the ground state and the excited states, can be well 
described by the LDA. This is an important character- 
istic to distinguish different phases driven by quantum 
fluctuations. Whereas we only test this approximation 
in ID, it is quite reasonable to conclude that it works 
also in 2D and 3D where fluctuations are less important. 
This generalization would greatly simplify many theo- 
retical studies on the actual superfluid since it builds a 
bridge connecting a confined system and a homogeneous 
one. This can also be used as a criterion to test the va- 
lidity of Bose-Hubbard model in the superfluid regime. 
In addition, the 3£7 resonance in the spectrum is a use- 
ful signature to characterize the density of the system. 
Therefore, the absorption spectroscopy is another exper- 
imental technique to study the density profile of the sys- 
tem besides the time-of-flight imaging. 

In conclusion, we have analysed in detail the dynami- 
cal response at zero temperature of the trapped ultracold 
bosons in an optical lattice subjected to lattice modula- 
tions. For the Mott-insulating system we identified sev- 
eral excitation processes. For the superfluid state, the 
presence of the harmonic trap induces slow- varying inho- 
mogeneities, which can be understood within the LDA. 
All these unique properties can be used examine whether 
the Bose-Hubbard model is a good realization of the ul- 
tracold atom system in a parabolic trapping potential. 
On the other hand, if one believes that the model can 
explain all the physics in this system, absorption spec- 
troscopy can be used as a powerful technique in revealing 
many basic features of the system, including its quantum 
state and density distribution. 
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